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Abstract 

Cell tracking algorithms which automate and systematise the analysis of time lapse im¬ 
age data sets of cells are an indispensable tool in the modelling and understanding of cellular 
phenomena. In this study we present a theoretical framework and an algorithm for whole cell 
tracking. Within this work we consider that “tracking” is equivalent to a dynamic reconstruc¬ 
tion of the whole cell data (morphologies) from static image datasets. The novelty of our work 
is that the tracking algorithm is driven by a model for the motion of the cell. This model may 
be regarded as a simplification of a recently developed physically meaningful model for cell 
motility. The resulting problem is the optimal control of a geometric evolution law and we dis¬ 
cuss the formulation and numerical approximation of the optimal control problem. The overall 
goal of this work is to design a framework for cell tracking within which the recovered data 
reflects the physics of the forward model. A number of numerical simulations are presented 
that illustrate the applicability of our approach. 

Keywords. Cell tracking, geometric evolution law, optimal control, phase field, finite elements. 


1 Introduction 

Cell migration is a fundamental process in cell biology and is tightly linked to many important 
physiological and pathological events such as the immune response, wound healing, tissue dif¬ 
ferentiation, metastasis, embryogenesis, inflammation and tumour invasion |jT|. Experimental 
advances provide techniques to observe migrating cells both in vivo and in vitro. Inferring dy¬ 
namic quantities from this static data is an important task that has many applications in biology 
and related fields. The field of cell tracking arose from this need and is concerned with the 
development of methods to track and analyse dynamic cell shape changes from a series of still 
images captured within a time frame (see for example for reviews). 

On the other hand, a major focus of current research is the derivation of mathematical mod¬ 
els for cell migration based on physical principles, e.g., |4|. Furthermore, such models appear 
to show good qualitative and quantitative agreement with experimental observations of migrat¬ 
ing cells. Despite this, very little research has focused on incorporating these mathematical 
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modelling advances into appropriate cell tracking algorithms. In a related work, we investig¬ 
ated fitting parameters in models for cell motility to experimental image data sets of migrating 
cells where observations of both the position of the cells and the concentrations of cell-resident 
proteins related to motility were available Q. 

In this study we present a first step towards the development of a framework for cell track¬ 
ing based on novel models of cell motility. Specifically, we propose a cell fracking algorifhm 
which can be fhoughf of as fiffing a simplified, yef physically meaningful, model for cell migra- 
fion fo experimenfal observafions and dafa. We focus on fhe selling, prevalenf in cell fracking 
problems, where only fhe posifion of fhe cell al a series of discrele limes is available and no 
furlher biological informalion is given. We presenl a malhemafical model based on physical 
principles for fhe cell movemenl fhal consisfs of a geomelric evolulion equafion. We Ihen for¬ 
mulaic an inverse problem, which fakes fhe form of a PDE conslrained opfimisalion problem, 
for fitting fhe model lo fhe sialic experimenfal observations. To solve fhe optimisation problem 
we propose an algorifhm based on fhe opfimal conlrol of geomelric evolulion laws |[^|7|. 

The objeclive of Ibis sludy is lo serve as a useful firsl sfep in fhe developmenl of cell frack¬ 
ing algorilhms in which fhe underlying model for fhe evolulion is based on physical principles, 
rafher fhan purely geomelric considerations. In fhis setting, one hopes lo attain eslimales of 
mofilily-relaled fealures such as Irajeclories, velocifies, persislence lenglhs, circularily, elc., 
which refiecl fhe physics underlying fhe model. We illuslrale fhe facl fhal fhe fracking pro¬ 
cedure we propose allows us lo incorporafe physically imporlanl aspecfs of cell migralion by 
including volume conservalion in fhe model for fhe evolulion. This is motivated by fhe observa- 
fion fhal, for many cells, while fhe surface area of fhe cell membrane may change significanlly 
during migralion Ihe volume enclosed by Ihe cell remains roughly conslanl Q. Of course 
olher physical aspecfs of fhe migralion could be included in fhe model, such as a sponlan- 
eous curvalure of Ihe membrane which is relevanl for more complex models of cell molilily 
involving Ihe Helfrich model Q. 

The remainder of our discussion proceeds as follows. In ^we briefly describe Ihe problem 
of cell fracking and inlroduce our approach lo cell fracking, which may be regarded as fitting a 
malhemafical model lo experimenfal image dafa sels. We presenl Ihe geometric evolution law 
model we seek to fit, which is a simplification of recently developed models in the literature 
that show good agreement with experiments | 4||8| - 131. We finish ^by reformulating our model 
into the phase field framework, which appears more suitable for the problem in hand, and we 
formulate the cell tracking problem as a PDE constrained optimisation problem. In ^ we 
propose an algorithm for the resolution of the PDE constrained optimisation problem and we 
discuss some practical aspects related to the implementation. In particular we note that the 
theoretical and computational framework may be applied directly to multi-cell image data sets 
and raw image data sets (of sufficient quality) without segmentation. In ^ we present some 
numerical examples for the case of 2d single and multi-cell image data sets. Einally in ^ we 
present some conclusions of our study and discuss future extensions and applications of the 
work. 


2 Problem Formulation 


2.1 Approaches to cell tracking 

Our focus is on developing whole cell tracking algorithms which track the morphology of the 
cell rather than particle tracking in which particles such as the cell centroid or cell resident 


proteins or (macro-)molecules are tracked 1141. A number of approaches have proved success¬ 


ful in cell tracking with level-set 115| or electrostatic based methods among the most widely 
used p^ . One feature of such methods is that the trajectories they generate are not physical 
in nature rather they are designed with the goal of achieving nice geometric properties, e.g., 
equidistribution of vertices, smoothness of the trajectories and so on. Our approach differs to 
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these purely geometric approaches in that we start with a model derived from physical prin¬ 
ciples and it is this model for the evolution that drives the tracking algorithm. In this sense our 
approach is similar in spirit to the parameter identification procedure described in Q as in both 
studies the goal may be regarded as fitting a mathematical model to experimental image data 
sets. 

We now summarise the main problems that must be addressed by a cell tracking algorithm. 
In general cell tracking consists of three main steps 

1. Segmentation: In this step the raw image data set is processed and the cells are separated 
from the background in each frame. 

2. Matching: The cells segmented in the first step must then be associated from frame to 
frame (note this is only relevant in the case of multiple cell image data sets) such that 
where possible (in practice cells may disappear or spontaneously appear in images) there 
is a one-to-one map that uniquely associates individual cells from one frame to the next. 

3. Linking: Finally the linking step consists of estimating dynamic data from the associated 
segmented static cells. 


In this work we will largely neglect the segmentation step. We assume either that we have 
segmented image data to work with or that the image data is of sufficient quality that the contrast 
between the cell and the background is clear and a simple thresholding step is sufficient to label 
the different cells. In the case of segmented image data, we assume this data consists of closed 
surfaces (or curves in 2d) that describe the boundaries of each individual cell. We briefly 
commenf on fhe role of segmenfafion in our approach and fhe facl fhaf if may be unnecessary fo 
exfracf fhe sharp contours of fhe cell boundaries in RemarkjT^ For ease of presenfafion we will 
also describe fhe algorifhm in fhe case of single cell image dafa and fhus fhe mafching sfep is 
redundanf. However fhe fheorefical aspecfs of our approach apply equally fo multiple cell dafa 
and one polenfially affracfive fealure of our cell fracking procedure is fhaf fhe mafching problem 
may be resolved implicifly wifhouf fhe need fo associafe cells form one frame fo fhe nexf, c.f., 
Remark 3.4 Our main focus is on fhe fracking sfep and in fhe remainder of fhis sfudy we oufline 
a fheorefical framework, a pracfical implemenfafion and numerical examples, for linking dafa 
af a series of discrefe times which allows fhe recovery of fhe whole cell morphologies in time. 


2.2 Model 

As menfioned above in confrasl fo many of fhe exisfing approaches for cell fracking, fhe frame¬ 
work we propose in fhis sfudy is based on tiffing a model, derived from physical principles, 
for fhe mofion of fhe cell fo experimenfal image dafa. The general class of models fo which 
our approach is applicable are PDF based models for fhe mofion, where fhe movemenf of fhe 
cell membrane is described by a geomefric evolufion law. We infroduce some nofafion for fhe 
formulafion of fhe model. 

We denofe by F fhe cell membrane, which is assumed fo be a closed smoofh orienfed d — 1 
dimensional hypersurface in d = 2,3, wifh oufward poinfing unif normal u. Given a 
funcfion t] defined in a neighbourhood of F, fhe fangenfial or surface gradienf of tj denofed by 
Vr is defined as 

Vrti := Vr/— Vt/• i/iv, (2.1) 

where V denofes fhe Cartesian gradienf in M'^. The Laplace-Belframi operator Ap is defined as 
fhe fangenfial divergence of fhe fangenfial gradienf, i.e.. 


Apr/ := Vr • (Vrr/). 

The mean curvafure H of T wifh respecf fo fhe normal u is defined as 

FT := Vr • 


( 2 . 2 ) 


(2.3) 
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In this study we model the evolution of the cell membrane as being governed by volume 
conserved mean curvature flow with forcing, given by 


V(x, t) = t) + r]{x, t) + Ay (f)) i/{x, t) 

r(o) = ro, 


on T{t),t G (0, T], 


(2.4) 


where T is the closed surface that represents the cell membrane, V is the material velocity 
of r, a is the surface tension and Ay(f) is a spatially uniform force accounting for volume 
conservation, physically this may be thought of as an interior pressure. The forcing function i] 
is the main driver of the directed migration. The model we present is phenomenological and 
hence it is difficult to directly relate r] to biophysical processes. However, as positive values of 
r] correspond to protrusive forces and negative values of r] correspond to contractile forces one 
interpretation of the forcing function rj is that it accounts for both protrusive forces generated 
by polymerisation of actin at the leading edge of the cell and contractile forces generated by 
the action of myosin motors at the rear of the cell. 


The evolution law (2.41 is a simplification of a large class of models that arise in the mod¬ 
elling of cell motility which take the following form 


V{x,t)= \^gi{H{x,t)) + g 2 {a{x,t)) + Xvit)ju{x,t) on r(f),f G (0,T], 

,r(o) = ro, 


(2.5) 


where gi models the dependence of the evolution on geometric quantities, such as resistance 


of the membrane to stretching which could be modelled by mean curvature terms as in or 
bending which can be modelled through the inclusion of Willmore or Helfrich flow type terms. 


The function g 2 appearing in (2.51 captures the dependence of the evolution on a vector of bulk 
and/or surface resident species a. The surface resident species a could, for example, satisfy 
another PDE such as a surface reaction-diffusion system 


d\^a + aVY'(^t)-V - D^Y{t)(^ = f {a) onr(f),iG (0,r], 
a(-, 0) = a°(-) on r(0), 


( 2 . 6 ) 


where a = (ai,..., a^^)^, is the number of chemical species involved, ai denotes the 
density of the ith chemical species, V is the material velocity of the surface. 


dl^a := dta + V ■ Va, 


(2.7) 


is the material derivative with respect to the velocity V, D is n diagonal matrix of positive dif¬ 


fusion coefficients and / (a) is the nonlinear reaction. Models of the form (2.51-(2.61 have been 


used successfully to model cell motility in Q ^ while models coupling evolution laws 


of the form (2.51 to bulk PDEs (i.e., equations posed in the cell interior) have been considered 


m 


Despite its simplicity the evolution law \2A) may be regarded as a prototype of the more 

. The geometric evolution component 


complex models for cell motility of the form (2.51 
( |2.5[ ) is often the most challenging component of the model to solve numerically and devel¬ 
oping an understanding of how to construct cell tracking algorithms assuming a geometric 
evolution law based model for the motion is an important first step towards developing track¬ 
ing algorithms based on more realistic physical models. In many applications it is also the 
case that the only information available from the data is the position of the cell membrane and 
no adequate model for the biochemistry of the motility related species involved is available. 
Without any knowledge of the relevant biochemistry it is difficult to identify which motility 
related species should influence the evolution let alone propose how the evolution depends on 


their distribution (i.e., a g 2 in (2.51) or a model for the species dynamics (i.e., an equation 


such as (2.6 1 ). Nevertheless one may still wish to extract dynamic quantities from static image 
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data sets, in this setting it may be reasonable to consider the evolution law (2.4 1 as a stand 
alone model for the motion as at least the mechanical aspects of the membrane evolution are 
accounted for through a physical model derived from basic physical principles. 


2.3 An optimal control approach to cell tracking 

The cell tracking approach we consider in this study corresponds to the following problem. 

2.4 Problem (Cell tracking). Given an initial cell membrane position and an observation 
of the position Tob^, find a space-time distributed forcing rj such that the evolution of the cell 
membrane, T{t),t € [0, T] satisfies (2.4 1 with r(0) = and r(T) the position of the cell 
membrane at time t = T, is close to Tofts. 

As the volume enclosed by the cell may vary over the images it is inappropriate to enforce 
conservation of a constant volume. Instead we enforce, with the help of a Lagrange multiplier 
Ay(f), that the volume enclosed by the cell is given by V{t) = + ^{Vobs — i-C- that 

the volume of the cell is a time-dependent linear interpolant of the volumes of the data. 


Problem 2.4 is an optimal control of a free boundary problem, where the free moving 
boundary problem is that of forced mean curvature flow and the control variable is the space 
time distributed forcing. The theory of optimal control of geometric evolution laws is in its 
infancy, in fact only recently has progress been made on the optimal control of parabolic equa¬ 
tions on evolving surfaces even in the case of prescribed evolution p^ . On the other hand 
the theory for the optimal control of semilinear parabolic equations is more mature (see, for 
example, (T^). We wish to exploit this fact and to this end we consider the phase field approx¬ 
imation of (2.41 given by the Allen-Cahn equation; 


Vip ■ un 


= Aip{x, t) - ^G'{ip{x, t)) - I {cGr]{x, t) - X{t)) in 0 x (0, T], 

= OonaO X (0,T], (2.8) 

= in O, 


where O C is a bulk domain, with normal uq, that contains T, is a diffuse interface 
representation of T*^ and e > 0 is a small parameter which governs the width of the diffuse 
interface. For details on the asymptotic analysis of (2 .81 and the convergence (as e —)■ 0) to 
a solution of p.4|) we refer the reader, for example, to ||20|-(2^ and references therein. The 


function G appearing in (2.8 1 is a double well potential, for example the quartic potential 


1 


G{p>) = - 1 )' 


(2.9) 


which has minima at ±1. The constant ^ G(r)^/^dr appearing in (2.8 1 is a scaling 


constant that depends on the double well potential. We enforce the time-dependent volume 


constraint following the approach of |211. Specifically our diffuse interface formulation of the 
constraint on the enclosed volume is given by a constraint on f^[ip(x, t)]_|_da;, where [a]+ = 
max(a,0). We define M^, the linear interpolant of f^[(p{x,t)]+dx of the initial and target 
diffuse interface data by 

M^{t) := j [v9°]+ -h ^ ([V’obsl-r - [<T°]+) d®, 


and determine X{t) in (2.8 1 such that M^p{t) = f^[ip{x, f)]+da;. We have used A (rather than 
Ay) for the Lagrange multiplier in (2.8 1 to reflect the fact that our constraint is on [p>{x, f)]+da;. 


However we shall refer to this constraint as a volume constraint in order to highlight the phys¬ 
ical feature the constraint is intended to model. We also investigated an alternative approach to 
enforcing the volume constraint via penalising deviations from a target volume following | [24] | 
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(see also Q), in our numerical tests this strategy proved less robust than the volume constraint 
proposed above. 

To formulate the cell tracking problem as a PDE constrained optimal control problem we 
define the objective functional we shall seek to minimise as follows 


= \ f i‘f{x,T)-<pobs{x)fdx + ^ [ [ r]{x,tfdxdt, ( 2 . 

^ JQ ^ Jo Jq 


10 ) 


where ipobs is a diffuse interface representation of the observation Tofts and 0 > 0 is a regular- 


isation parameter. The first term on the right of ( |2.10[ ) is the so called fidelity term that measures 
the distance between the solution to the model and the target data and the second term is the 


regularisation which is necessary to ensure a well-posed problem (for example see 1191). 

Our optimal control approach to the cell tracking problem may now be stated as the follow¬ 
ing minimisation problem. 


2.5 Problem (Optimal control problem). Given an initial diffuse interface representation of the 
cell membrane position and an observation of the position (fobs-: find a space-time distributed 
forcing tf : O x [0, T] —)• M such that with a solution of (^1 with initial condition ip{-, 0) = 
the forcing rj* solves the minimisation problem 


min J((p, rj), with J given by p.lOj). 
n - 


( 2 . 11 ) 


2.6 Optimality conditions 


To apply the theory of optimal control of semilinear PDEs for the solution of the tracking prob¬ 
lem, we briefly oufline fhe derivafion of fhe opfimalify condifions, for furfher defails see for 
example | fT9l[25| . Infroducing fhe Eagrange mulfiplier (adjoinf sfafe) p, we define fhe Eag- 
rangian functional 


[ [ (dtip{x,t)-Aip{x,t) 

Jo Jn \ 


+ ^G'{(f{x,t)) + ^{cG'n{x,t) - \{t)) jp(a;,f)d®df. 


( 2 . 12 ) 


Requiring sfafionarify of fhe Eagrangian wifh respecf fo fhe adjoinf sfafe yields fhe sfafe equa- 


fion (2.8 1 and requiring sfafionarify of fhe Eagrangian, af fhe opfimal confrol tj* and associafed 
optimal sfafe ip*, wifh respecf fo fhe sfafe and fhe confrol, yields fhe (formal) firsf order opfim¬ 
alify condifions 


6^C{p*,r]*,p)p = 0, \/p : p{x,0) = 0, (2.13) 

6rjC{p*,r]*,p)rj = 0 , V??. (2.14) 


Condition ( 2.13| ) yields 
fhe adjoinf sfafe p, 


fhe adjoinf equafion, which is fhe following linear parabolic PDE for 


dtp{x,t) =-Ap{x,t) + ^G"{p{x,t))p{x,t) inDx(0,r], 

< Vp-i^n =0 onc)Dx(0,T], (2.15) 

^p{x,T) = p{x,T) - pobs{x) inD. 


Nofe fhaf equation ( |2.15| ) is posed backwards in fime and hence is equipped wifh ferminal con¬ 
difions. Condifion (2. 14|) fogefher wifh fhe Riesz represenfafion fheorem yields fhe opfimalify 


condition (c.f., 1191) 


Sr,C{p*,Tl*,p) = 07]* + ^p = 0 . 


(2.16) 
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2.7 Remark (Choice of the potential). We note that our approach to the optimal control prob¬ 
lem involving the formulation of the adjoint problem appears to require a smooth potential G 
(c.f., p.9[ )). The formulation of the adjoint problem is to our best knowledge an open prob¬ 
lem for other widely used, but non smooth or unbounded, potentials such as the obstacle or 
logarithmic potential. 


3 Practical considerations, implementation and algorithm 


As is standard, we use the optimality conditions to construct an iterative optimisation loop to 
solve the optimal control problem, Problem [2.5| The basic idea is that in each step of the loop 
we first solve the state equation ( |2.8| ) with a given control, then solve the adjoint equation ( |2.15| ) 
with the computed states and then update the control using the optimality condition \2.\6\ . For 
this initial study to ensure robustness of the algorithm and to aid in the clarity of the exposition 
we employ a simple gradient based update of the control | |2^ . Given rj^ and p we compute 
the updated control via steepest descent. That is we choose as an update direction the 
negative gradient, the formula for the update of the control is 


+ ^p{x,t)j , (a;, f) G 12 x [0, T), (3.1) 


where a is a step size. For simplicity in this study we take a constant step size of a = 0.01. 

For the termination criteria for the algorithm we stop if the objective functional J is less 
than a given tolerance tolj, the update in the control is less than a given tolerance, i.e., if 

or if a maximum number of iterations Kmax is reached. In 


practice the forward ( 2.81 and adjoint equations (2.151 must be approximated and to do this we 
utilise a finite element method, using continuous piecewise linear elements. The details of the 
numerical method for the approximation of the forward and adjoint equations are given in[A| 
The cell tracking algorithm we propose may now be stated in pseudocode as follows (the 
notation employed is as in[A]); 


3.1 Optimal control based cell tracking algorithm 

Require: Data: {p^obs)h the initial and target (discrete) diffuse interface data. 

Numerical Parameters: T > 0 end-time and M > 0 number of timesteps. 

Optimisation Parameters: tolerances tolj, tolrj, and Kmax- 
Initial guess for the control: {rjh)^ ■= iVh)^ G V, z = 0,..., M. 
k := 0 

while ( a{9{rihf + > tolr^, J > tolj and k < Kmax) do 

V e n ' L2(nx[0,T)) ' / 

Solve state equation for (A*)^"''^}, z = 1,..., M with and initial data 

= >c.f.,0 


Solve adjoint equation for = M — 1,..., 0 with computed and with 

terminal data - {^obs)h- 

Update control = {r]D^ — a{9{p\)^ -|- z = 0,... ,M. 


Compute J according to (2.101. 

k := k + 1. 


end while 
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3.2 Remark (Segmentation and image data). An important aspect of any cell tracking al¬ 
gorithm is its ability to extract data suitable for the tracking algorithm from the experimental 
image data set. In many cases the experimental image data set is grayscale data with the in¬ 
tensity (brightness) indicating whether a point is in the interior of a cell, i.e., points inside the 
cell appear bright for example and points outside appear dark. For many tracking algorithms 
this intensity data is then post processed via a segmentation algorithm (e.g., active contour 
methods |27 ^) to yield sharp interface representations of the cell membrane. Assuming a 
sharp interface representation of the cell membrane is available, diffuse interface representa¬ 
tions may be easily initialised (see for example Q). We note however that the raw intensity 
data produced by many imaging procedures may already be close to a diffuse interface rep¬ 
resentation of the cell, this is typically the case when the data is relatively free of noise and 
the contrast between the cell and the background is high. In this case one may wish to exploit 
this fact in the algorithm and work with the raw image data set itself (or a post processed e.g., 
thresholded version), thus circumventing the extra error induced by segmentation. 

3.3 Remark (Observations at multiple points in time). For clarity of exposition we focus on the 
case of fitting to a single observation. The approach generalises straightforwardly to multiple 
observations with the first term in ( 2.10[ ) simply replaced by a sum over the distinct times at 
which the observations are taken of the difference between the solution (at the appropriate time) 
and the target data. 


3.4 Remark (Multiple cells and matching problems). As mentioned above a major focus of 
many cell tracking algorithms is to track multiple cells in the same image and the resolution of 
the so called matching problem. Our approach can be applied to multi-cell image data. Here 
and t^ohs would be diffuse interface representations of the multi-cell image data set and 
the diffuse interfaces would consist of multiple disjoint phases. The remaining aspects of the 
approach remain unchanged and the matching problem is solved implicitly in the computation 
of the optimal control. 

There are however multiple practical issues which arise in this setting related to the separa¬ 


tion between distinct cells, which affects the choice of e, and the fact that the evolution law ( 2.81 


allows changes in the topology of the phases which may lead to cell splitting, the annihilation 
of a phase (which would correspond to the disappearance of a cell) or the nucleation of a phase 
(i.e., the spontaneous appearance of a cell). We intend to comment on practical approaches to 
multi-cell tracking elsewhere. 


4 Numerical examples 


We now present some benchmark numerical examples illustrating the application of the al¬ 
gorithm to artificial image data sets. For all the simulations we report on in this section, in 


the state equation (2.8 1 we set e = 0.1, and we took the end-time T = 0.4. As mentioned 
previously the parameter e governs the width of the diffuse interface and should in general 
be taken as small as is computationally feasible, smaller values of e necessitate a finer grid, 
in this initial study with uniform grids we set e = 0.1 as the CPU times become prohibitive 
for smaller values of e. The end time T corresponds to the nondimensionalised time between 
snapshots and could in principle be related to an acquisition time between images given real 
biological data. For each of the experiments, apart from those of ^4.4[ we set the initial guess 
for the control to be constant in space and time (zero in the single cell case and one for the 
multi-cell examples). For the approximation of the forward and adjoint equations we used 
triangulations with 8321 DOFs in all the simulations, apart from those of ^4.3[ and selected a 
uniform timestep r = 1 x 10“^. The same numerical parameters for the optimisation algorithm 
were used for all the experiments and are given in Table [T] In every example we report on the 
algorithm terminated due to the update of the control being less than the prescribed tolerance. 
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a 9 

tolj 

toljj 

K 

0.01 0.01 

1 X 10“^ 

1 X 10-4 

3500 


Table 1: Parameter values used for the numerical simulations. 

The technical details of the hardware used to carry out the simulations is given in Remark 

oi 

4,1 Remark (Hardware details). All the numerical experiments have been performed on the 
high performance cluster (HPC) at the University of Sussex. Each of the simulations was 
carried out in serial using a single core of the cluster. The HPC cluster currently consists of 
3140 cores with an even mixture of Intel and AMD CPUs. The majority of the cluster are 64 
core AMD nodes with 256GB RAM per node, and a smaller number of 512GB RAM nodes. 
The cluster uses the high-performance Lustre clustered-filesystem for EO, and currently stands 
at 298TB of storage for research use. 


4.2 Application to synthetic data 


Here we apply the algorithm to a single synthetic cell data set taken from the PhagoSight web¬ 
site http : //WWW. phagosight. org/synData . php. The synthetic cell was generated 
as a mixture of Gaussians with Poisson noise that varied over time to simulate the displacement 
and change of shape of a neutrophil as observed in a Zebrafish embryo. The data for analysis 
consisted of points on the synthetic cell membrane at a series of times (for simplicity we used 
2d data, i.e., the cell membrane was a Id curve embedded in M^). The initial and target curves 
we took as test data for the algorithm are shown in Figure [T(a)| To apply our algorithm, based 
on diffuse interface representations, we define fhe domain U:=[0, 8] x [0,6] which was such fhaf 
bofh fhe initial and fargef curves were confained in fhe domain. We fhen consfrucfed diffuse in- 
ferface represenfafions of fhe fargef dafa following fhe procedure described in Q. Figures[T(b)| 


and 1(c) show fhe diffuse inferface represenfafions of fhe initial and fargef dafa respecfively. In 
order fo invesfigafe fhe influence of fhe volume consfrainf on fhe compufed cell morphologies 
we performed fwo experimenfs, in fhe firsl we simply considered fhe forced Allen-Cahn model 
for fhe evolution wifh no volume consfrainf, i.e., ( |2.8| ) wifh A = 0, and in fhe second we in¬ 
cluded fhe volume consfrainf as described in ^ The algorifhm fook 1996 iferafions fo meef 
fhe stopping criferia wifh no volume consfrainf and 2479 iferafions wifh fhe volume consfrainf, 
corresponding fo CPU times of 25238 and 119216 seconds respecfively. 

Figureshows fhe value of fhe objecfive funcfional againsf fhe number of iferafions of fhe 
opfimisafion algorifhm wifh and wifhouf fhe volume consfrainf. We observe similar behaviour 
in bofh cases wifh an initial rapid decrease in fhe objecfive funcfional followed by a more 
gradual reduction wifh each iferafion as we approach fhe minimum. Figure shows fhe zero 
level-sef of fhe compufed solufion using fhe opfimal confrol af fhe final fime wifh and wifhouf 
fhe volume consfrainf. The curve corresponding fo fhe zero level sef is shaded by fhe value of 
fhe compufed opfimal confrol. The background shading corresponds fo fhe fargef dafa. In bofh 
cases fhe posifion of fhe zero level-sef of fhe compufed solufion shows good agreemenf wifh 
fhe fargef dafa. Qualifafively we observe cells wifh a clearly defined “fronf” and “rear”, wifh 
fhe compufed confrol corresponding fo profrusive forces af fhe fronf and confracfive forces af 
fhe rear. 

Figure shows fhe area enclosed by fhe zero level-sef of fhe solufion wifh fhe opfimal 
confrol wifh and wifhouf fhe volume consfrainf fogefher wifh fhe linear inferpolanf of fhe areas 
of fhe dafa. We see fhaf wifhouf fhe volume consfrainf fhe area initially decreases fhen rapidly 
increases as we approach fhe final fime whilsf wifh fhe volume consfrainf (nofe fhe consfrainf 
is acfually on fhe mass rafher fhan fhe volume) fhe area is close fo fhe linear inferpolanf of fhe 
areas of fhe dafa. In ferms of fhe compufed cell morphologies. Figure j^shows snapshofs of fhe 
compufed cell membranes (zero level-sefs) for fhe fwo differenf cases. We clearly observe fhaf 
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(a) Initial (red curve) and target (green curve) synthetic data. The cell 
centroids are shown together with the trajectory of the linear interpolant 
of the cell centroids (black line). 



0.0 10 2.0 3.0 4.0 5.0 6.0 7.0 8.0 0.0 10 2.0 3.0 4.0 5.0 6.0 7.0 8.0 

(b) Initial data (c) Target data ((fobs)- 


Figure 1: Initial and target data for the example of §4.2 


the intermediate snapshot (blue curve) encloses a much smaller area if the volume constraint is 
not included in the algorithm. In Figurewe report on the trajectory and speed (magnitude of 
the velocity) of the centroid (center of mass) of the zero level-set of the computed solution with 
the optimal control, with and without the volume constraint. We observe similar trajectories 
with and without the volume constraint and in both cases we observe an increase in the speed 
as we approach the end-time. However in the case of no volume constraint this increase is 
more marked with a sharp spike in the centroid velocity observed close to the final time. The 
increasing centroid speed we observe may be unphysical and if a (roughly) constant centroid 
velocity is desired one strategy may be to impose pointwise constraints on the control, this 
would prevent the large increase in the maximum and minimum values of the control observed 
during the simulations as we approach the final time as shown in Figure Another possible 
strategy would be to modify the regularisation in ( |2.10| ). 

Figurej^shows the fidelity term ||(/9obs(®) — ^)^||l 2 ( 0 ) without the volume 

constraint versus the number of iterations (where k corresponds to the optimisation iteration 
number). The fidelity term may be considered as a quantitative measure for the “goodness of 
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(a) Without the volume constraint (b) With the volume constraint 

Figure 2: The value of the cost functional versus the number of iterations for the experiments of 
§ |4.2| with and without the volume constraint. We observe a rapid decrease in the cost initially 
followed by a much more gradual decrease as we approach the minimum, this is as expected since 
the steepest descent algorithm is used for the update of the control. 
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(a) Without the volume constraint 
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(b) With the volume constraint 


Figure 3: Zero level-set of the solutions {(p{x,T)) computed using the approximated optimal 
control (r]*) with and without the volume constraint for the experiments of |4.2[ The curve (zero 
level-set of T)) is shaded by the approximated optimal control T)) and the background 
by the target data (cpobsix)). The color-bar corresponds to the scale for r]*{x,T). We see good 
agreement between the zero level-set of the data computed with the optimal control and the target 
data in both cases. 

fit" of the computed data to the observations. We observe a steady decay in the fidelity term as 
we approach the optimal control in both cases. 

4.3 The effect of mesh refinement 

In this section we investigate the effect of the mesh-size on the results by refining the mesh 
whilst keeping the time step r constant. We report on the value of the fidelity term computed 
using 0*, i.e., the forward state computed with the optimal control. The initial and target data 
are taken to be the same as in ^4.2| and we employ the algorithm with the volume constraint. 

The initial guess for the control is taken to be zero in each case. The results from the mesh 
refinements are presented in Table We observe a reduction in the fidelity term as we refine 
the mesh which implies an improved fit to the observed data. Although in principle it would 
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Figure 4: Area enelosed by the eell for the experiments of §4.2 with and without the volume 


eonstraint. The eell shrinks eonsiderably during the evolution without the volume eonstraint whilst 
a good fit to the linear interpolant of the area enelosed by the data is observed with the volume 
eonstraint. 




Figure 5: Zero level-sets of the solutions eomputed t)) with the optimal eontrol {ri*{x, t)) for 
the experiments of ^4. 2 [ with and without the volume eonstraint at f = 0 (red), t = 0.35 (blue) and 
f = 0.4 (green). We observe that the volume enelosed by the blue eurve is signifieantly smaller 
than the volumes enclosed by the red and green curves without the volume constraint whilst this is 
not observed if the volume constraint is included. 

be interesting to investigate the influence of refining both the timestep and mesh-size on the 
computed results, our tests indicate that the algorithm breaks down for time steps significantly 
larger than 0.01 and hence, refinement of the timestep and mesh-size together becomes com¬ 
putationally prohibitive. 
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(a) Centroid trajectories together with centroids of the data. 



Figure 6: Trajectories and speeds of the centroid of the zero level-sets of the solution with the 
optimal control with and without the volume constraint for the example of |4.2 


DOFs 

\\^{x,t) - ^obs\\L 2 {n) 

CPU time (secs) 

545 

0.348165 

3142 

2113 

0.176357 

26050 

8321 

0.154305 

119216 

33025 

0.135178 

291597 


Table 2: Mesh refinement. 


4.4 The influence of the initial guess for the control 

Here we apply the algorithm with the volume constraint on the simple example of a translated 
circle to illustrate the effect that the choice of the initial guess for the control rj has on the 
solution of the problem. To apply our algorithm we define the domain 17 to be [—3,6] x [—3,3] 
with a trianguiation of 8321 grid points. We seiected a uniform timestep r = 1 x 10“^ and 
set the interfaciai thickness e = 0.1. We took the end-time T = 0.8. The remaining numerical 
parameters for the optimisation algorithm are as given in Table [T] The initial data was taken to 
be a smoothed (by running a few steps of the Allen-Cahn solver) version of the function taking 
the value 1 inside Hi (0,0) (a circle of radius 1 centred at the origin) and -1 in 17/Hi(0,0). 
The target data was taken to be a smoothed (by running a few steps of the Allen-Cahn solver) 
version of the function taking the value 1 inside Hi (3,0) and -1 in 17/Hi(0, 0). Figurej^shows 
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(a) Minimum value of the control for a series of iterations against 
time for the algorithm without volume constraints. 



(b) Maximum value of the control for a series of iterations against 
time for the algorithm without volume constraints. 



(c) Minimum value of the control for a series of iterations against 
time for the algorithm with volume constraints. 



(d) Maximum value of the control for a series of iterations against 
time for the algorithm with volume constraints. 


Figure 
for the 


7: Minimum and maximum values of the control rj with and without the volume constraint 


example of §4.2 
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(a) Without the volume constraint 
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(b) With the volume constraint 


Figure 8: The value of the fidelity term versus the number of iterations for the experiments of §4.2 
with and without the volume constraint. 


the initial and target diffuse interface data. To illustrate the effect of the choice of initial guess 



(a) Initial data ((^°). 


(b) Target data {(fobs)- 


Figure 9: Initial and target data for the examples of §4.4 


on the algorithm, we consider two different values for the initial guess, firstly we set r/ = 0 
and secondly we set r/ = c • V(/?, where c = (2.5,0), i.e., in the latter case the initial guess 
depends on the solution to the Allen-Cahn equation. In both cases we used the algorithm with 
the volume constraints. With the zero initial guess the algorithm took 3262 iterations to meet 
the stopping criteria corresponding to a CPU time of 320433 seconds. With the second choice 
of initial guess the algorithm took 2056 iterations to meet the stopping criteria corresponding 
to a CPU time of 228173 seconds respectively. 

Figure shows the zero level-set of the computed solution using the optimal control at the 
final time. The curve corresponding to the zero level-set is shaded by the value of the control 
with the background shading corresponding to the target data. In both cases the position of 
the computed curve (zero level-set) with the optimal control shows good agreement with the 
target data. Figure shows snapshots of the computed zero level-sets with the two different 
initial guesses. For the case with the initial value of ry = 0, we observe in Figure 11(a) that the 
interface remains close to the initial position for most of the time of the simulation, and at the 
very last moment it shrinks to a point with a new phase nucleated at the position of the target 
data corresponding to a change in topology. With the second choice of initial guess (ry = c-Vip) 


we observe in Figure 11(b) that there is a gradual motion towards the target position with no 
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changes in topology. Figure [T^ shows the area enclosed by the zero level-set of the computed 
solution with the optimal control with the two different initial guesses together with the linear 
interpolant of the areas of the data. We observe a sharp increase in area towards the end of the 
time interval with the zero initial guess as the new phase is nucleated. With the second choice 
of initial guess, the area of the computed curve exhibits a good fit to the linear interpolant of 
the areas of the data. 



(a) With initial guess 77 = 0 


(b) With initial guess rj = c-Vip 


Figure 10: Zero level-set of the solutions {ip{x,T)) eomputed using the approximated optimal 
eontrol t)) for the experiments of ^4.4 The eurve (zero level-set of (p{x, T)) is shaded by 
the approximated optimal eontrol {r]*{x, T)) and the baekground by the target data {ipobs{x)). The 
eolor-bar eorresponds to the seale for 7f{x, T). We see good agreement between the zero level-set 
of the data eomputed with the optimal eontrol and the target data in both oases. 



(a) With initial guess rj = 0 


(b) With initial guess rj = c-Vip 


Figure 11: Zero level-sets of the solutions eomputed ((p{x, t)) with the optimal eontrol (rj* {x,t)) for 
the experiments of § |4.4| at t = 0 (red), t = 0.2 (blaok), t = 0.6 (blue), t = 0.7 (orange), t = 0.789 
(pink) and t = 0.8 (green). We observe the nuoleation of a phase and a ohange in topology with the 
zero initial guess whilst there are no evident ehanges in topology and the zero level-set maintains a 
fixed topology in the ease of the nonzero initial guess. 
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Figure 12: Area enclosed by the curve for the experiments of ^4.4 A good fit to the linear 


interpolant of the areas is only observed with the nonzero initial guess. We observe a rapid increase 
in the area near the end time for the zero initial guess, this corresponds to the time at which a new 
phase is nucleated, c.f.. Figure [TT(a)j 


4.5 Application to multi-cell image data sets 

We now apply the algorithm to the case of multi-cell image data sets. As a proof-of-concept 
we consider the simplest possible scenario where we have an initial and desired dataset both 
consisting of two cells that are well separated. 

For the first experiment we defined the initial data and target data as follows. Defining the 
domain D to be [— 2 , 8 ] x [— 2 , 2 ] we defined the subdomains Di, ^ 2 ) ^3 and D 4 to be the 
simply connected bounded domains with boundary curves Fi, r 2 , Fa and F 4 defined by (the 
curves Fi, F 2 and F 3 and F 4 are the zero level-sets of the diffuse interfaces shown in Figures 
13(a)| and |13(b)| respectively). 


Fi := |x G D j Xa -|- *2 — 0.8^ -|- 0.1 sin(4xi) -|- 0.1 sin(3x2) = 0} , 

F2 := |x G D 1 - 2^ + {x 2 - 0.6)^ - 0.7^ -1-0.1 sin + 0-3sin(2x2) = > 

F3 := jx G D j (xi — 0.4)^ -|- (x2 — 0.5)^ — 0.8^ -|- 0.1 sin(6xi) -|- 0.1 sm(7x2) = O} , 

F4 := |x G Vt \ — 2.5^ -|- (x2 — 1)^ — 0.7^ -|- 0.1 sin sin(1.5x2) = o| . 

We then set the initial and target data to be a smoothed (by running a few steps of the Allen- 
Cahn solver) version of the function 


7 ^° = 


1 for X G U O2, 

—1 for X G D/ (Di U ^2) 


and 


^obs 


1 for x G D3 U D4, 

— 1 for a; G D/ (fls U D4). 


Figure [T^ shows the initial and target diffuse interface data. 

As previously, we compare the results of the algorithm with and without the volume con¬ 
straint. For this experiment, the algorithm took 2035 iterations to meet the stopping criteria 
with no volume constraint and 2199 iterations with the volume constraint, corresponding to 
CPU times of 28608 and 105750 seconds respectively. 
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Cost 


Figure [T^ shows the value of the objective functional against the number of iterations of 
the optimisation algorithm with and without the volume constraint. Figure shows the zero 
level-set of the computed solution using the optimal control at the final time with and without 
the volume constraint shaded by the value of the control with the background shading corres¬ 
ponding to the target data. The results are similar to the single cell simulations of §4.2| with 
an initial rapid decrease in the cost followed by a subsequent gradual decrease. The cells (zero 
level-sets) computed with the optimal control show good agreement with the target data for both 
versions of the algorithm and for both cells. For each of the versions of the algorithm, both of 
the computed cells again posses a clearly defined “fronf” and “rear” similar fo fhe single cell 
case. 

Figure [T^ shows fhe area enclosed by fhe zero-level sef of fhe compufed solufion wifh fhe 
optimal confrol wifh and wifhouf fhe volume consfrainf fogefher wifh fhe linear inferpolanf 
of fhe areas of fhe dafa. We observe analogous behaviour fo fhe single cell. In terms of fhe 
compufed cell morphologies, Figure [T7] shows snapshofs of fhe computed zero level-sefs for 
fhe fwo differenl versions of fhe algorifhm. We see fhaf in fhis multi-cell selling fhe algorifhm 
has implicilly solved fhe mafching problem by generaling fwo disjoin! cells whose lopology 
remains fixed fhroughoul fhe evolulion. We observe fhaf fhe loss of volume in fhe case of no 
volume consfrainf corresponds fo one of fhe cells in fhe inlermediale snapshol (blue curve) 
enclosing a much smaller area. 



Figure 13: Initial and target data for the examples of §4.5 



(a) Without the volume constraint (b) With the volume constraint 


Figure 14: Cost functional versus the number of iterations for the examples of §4.5 


4.6 An example with topological change 


Of course, in general our algorithm may generate cells whose topology is not fixed as in §4.4 
To this end we report on another experiment. 
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(a) Without the volume constraint 


(b) With the volume constraint 


Figure 15: Zero level-set of the solutions {ip{x,T)) eomputed using the approximated optimal 
eontrol {rj* {x,t)) with and without the volume eonstraint for the examples of |4.5 The eurve (zero 
level-set of (p{x, T)) is shaded by the approximated optimal eontrol {r]*{x, T)) and the baekground 
by the target data ((pobsix)). The eolor-bar eorresponds to the seale for r]*{x,T). We see good 
agreement between the zero level-set of the data eomputed with the optimal eontrol and the target 
data in both oases. 



Figure 16: Area enclosed by the cell for the experiments of ^4.5 with and without the volume 


constraint. As with the single cell data, the area (now the sum of the area of the two cells) shrinks 
considerably during the evolution without the volume constraint whilst a good fit to the linear 
interpolant of the area enclosed by the data is observed with the volume constraint. 


Defining the domain D to be [—2, 6.3] x [—2.5, 2.5] we defined the subdomains Di, D 2 , D 3 
and D 4 to be the simply connected bounded domains with boundary curves Fi, r 2 , Fs and F 4 
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(a) Without the volume constraint (b) With the volume constraint 


Figure 17: Zero level-sets of the solutions computed with the optimal control for the examples of 
§ |4.5| with and without the volume constraint at f = 0 (red), t = 0.35 (blue) and t = 0.4 (green). 
The volume enclosed by both cells shrinks during the evolution without the volume constraint 
whilst this is not observed if the volume constraint is included. Both with and without the volume 
constraint, the implicit solution of the matching problem in this case generates two disjoint cells 
which do not change in topology. 


defined by (see Figure [T8]) 

Fi := \^x ^ Q \ xf + X 2 — 0.9^ -|- 0.1 sin(4.5xi) -|- 0.11 sin(3x2)) = O} , 
r2 := G n I (xi — 5)^ + X 2 — 0.7^ -1-0.1 sin sin(2x2) ~ ’ 

Tg := {x G 0 I (xi - 0.35)2 q _ q g 2 q.I sin( 6 xi) -h 0.1 sin( 7 x 2 ) = O} , 

, /_ 1 i\2 n ^2 r, 1 


r4 := G n I (xi — 0.3) -|- (x 2 — l.ly — 0.72 — O.l sin + O-l sin(l.5x2) ~ ^ j • 

We then set the initial and target data to be a smoothed (by running a few steps of the Allen- 
Cahn solver) version of the function 


Q j I for * G Hi U H2, 

^ [ —I for a; G H/ (Hi U H2) 


and 


^obs — 


I for a? E ils U ^^4, 

— I for a; G H/ (H 3 U H 4 ). 


Figure [T^ shows the initial and target diffuse interface data. 

As previously, we compare the results of the algorithm with and without the volume con¬ 
straint. For this experiment, the algorithm took I960 iterations to meet the stopping criteria 
with no volume constraint and 1937 iterations with the volume constraint, corresponding to 
CPU times of 27553 and 93150 seconds respectively. 

Figure [T^ shows the value of the objective functional against the number of iterations of 
the optimisation algorithm with and without the volume constraint. Figure shows the zero 
level-set of the computed solution using the optimal control at the final time with and without 
the volume constraint shaded by the value of the control with the background shading corres¬ 
ponding to the target data. The results are similar to the previous simulations with an initial 
rapid decrease in the cost followed by a subsequent gradual decrease and good agreement with 
the target data for both versions of the algorithm and for both cells. For each of the versions of 
the algorithm, both of the computed cells again posses a clearly defined “front” and “rear”. 

Figure shows the area of the domain in which the computed solution is positive with 
and without the volume constraint together with the linear interpolant of the areas of the data. 
For this experiment we observe that the area enclosed by the computed cells differs signific¬ 
antly from the linear interpolant areas of the data both with and without the volume constraint. 
This may be due to the change in topology of the interface during the evolution. Figure 
shows snapshots of the computed zero level-sets for the two different versions of the algorithm. 
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Unlike the previous examples we see that for this particular choice of initial and target data, 
the algorithm yields cells which change in topology with one of the curves shrinking until it 
disappears whilst the other splits eventually becoming two disjoint curves. Thus our algorithm 
generates trajectories corresponding to the annihilation (via shrinking) of one cell whilst the 
other cell splits to form the two cells observed in the image data set. 



Figure 18: Initial and target data for the examples of §4.6 



(a) Without the volume constraint (b) With the volume constraint 


Figure 19: The value of the eost funetional versus the number of iterations for the examples of 
§ |4.6| with and without the volume eonstraint. We observe a rapid deerease in the eost initially 
followed by a mueh more gradual deerease as we approaeh the minimum, this is as expeeted sinee 
the steepest deseent algorithm is used for the update of the eontrol. 


4.7 Comments on the numerical experiments 

The CPU times for each of the experiments is of the order of hours. For all the experiments 
the number of iterations required before the stopping criteria is met are similar, however this 
leads to simulations with the volume constraint taking approximately four times as long (in 
terms of CPU time) as those without the volume constraint. This is due to the iterative nature 
of the algorithm used to compute the Lagrange multiplier c.f., which necessitates multiple 
solves per timestep. We note that the CPU times of the algorithms may be too large for many 
applications. In light of this we mention that the stopping criteria we have used may be too 
strict for many applications and that a significant decrease in the cost function together with 
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(a) Without the volume constraint 


(b) With the volume constraint 


Figure 20: Zero level-set of the solutions {^p{x,T)) computed using the approximated optimal 
control (t]* {x,t)) with and without the volume constraint for the examples of |4.6 The curve (zero 
level-set of ip{x, T)) is shaded by the approximated optimal control ir]*{x, T)) and the background 
by the target data (^Pobsix)). The color-bar corresponds to the scale for t]*{x,T). We see good 
agreement between the zero level-set of the data computed with the optimal control and the target 
data in both cases. 


a reasonable fit to the data is observed after as few as 50 iterations (which reduces the CPU 
time by a factor of around 40) and that for many applications this level of fit may be sufficient 
hence the stopping criteria could be relaxed. Finally we mention that the current solution pro¬ 
cedure based on uniform grids and serial solution of the forward and adjoint problems may be 
improved and we are currently investigating combining adaptive grids with a parallel solver for 
the forward and adjoint problems which gives a significant speed up but presents new technical 
challenges which we wish to avoid in this paper to maintain clarity of exposition. 


5 Conclusion 

In this study we presented a first step towards the development of cell tracking algorithms 
based on physical models for cell migration. The presented algorithm seeks to track whole cell 
morphologies and is applicable to single cell or multi-cell image data sets. Our approach may 
be regarded as a model fitting procedure in which a physically derived model for the evolution 
of the cell or cells is fitted to experimental image data sets. The algorithm is based on the 
theory of optimal control of PDFs and full details of the derivation and implementation of 
the algorithm are given. We also present a number of numerical experiments illustrating the 
performance of the algorithm with synthetic representative single cell and multi-cell image data 
sets. 

The key novelty of our approach is that the model for the evolution of the cell (or cells), 
which drives the tracking procedure, is based on a relevant simplification of existing physically 
derived models for cell motility that reproduce many experimentally observed aspects of cell 
migration (e.g., |4|). Thus this study is an important step towards the development of cell track¬ 
ing algorithms in which the recovered trajectories are physically meaningful. This is in contrast 
to the majority of existing algorithms for whole cell tracking in which the models for the evol- 
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Figure 21: Area enelosed by the eell for the experiments of §4.6 with and without the volume 


eonstraint. As with the single eell data, the area (now the sum of the area of the two eells) shrinks 
eonsiderably during the evolution without the volume eonstraint, whilst the ineorporation of the 
volume eonstraint yields a better fit to the linear interpolant of the areas. Unlike the previous 
examples however, even with the volume eonstraint the fit to the linear interpolant of the areas of 
the data is poor. 





(a) Without the volume constraint 


(b) With the volume constraint 


Figure 22: Zero level-sets of the solutions eomputed with the optimal control for the examples of 
§ |4.6| with and without the volume constraint at f = 0 (red), t = 0.3 (blue), t = 0.396 (black) 
and t = 0.4 (green). We observe that in this case the difference between the two schemes is less 
pronounced. It is also clear that both with and without the volume constraint, the implicit solution 
of the matching problem in our algorithm in this case leads to the annihilation of one cell (as it 
shrinks to a point) while the other cell splits with the zero level-set changing in topology from a 
single closed curve to two disjoint closed curves. 


ution of the cell that underly the tracking procedure are purely geometric in nature neglecting 
completely the physics of cell migration | T5|[T6p9 |. One significant advantage of the approach 
to cell tracking we propose, is that the physics of the model driving the evolution of the cell are 
reflected in the recovered dynamic data. Thus it is possible to encode physical features of cell 
migration into the tracking procedure. We illustrate this fact by including volume conservation 
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in the model. Comparing the results of the tracking algorithm with and without volume con¬ 
servation, we observe, that in a number of simulations neglecting volume conservation leads 
to physically unrealistic cell morphologies with a significant reduction of the cell volume in 
the recovered morphologies whilst this undesirable effect is no longer evident if volume con¬ 
servation is included. We note that volume conservation is more physically relevant in three 
dimensions. This is since large changes in volume in two-dimensional imaging data of cells 
(i.e., the area of the projections of the cell on to a two-dimensional plane) are often observed in 
experimental data despite the cells conserving their enclosed (three-dimensional) volume. 

Of course volume conservation is only an example of the kind of biology or physics one 
may wish to encode in the algorithm. A number of models that fit into our framework have 
been proposed incorporating more complex biophysics such as spontaneous curvatures Q, for 
Helfrich type models, adhesion for the migration of cells on substrates or in the ECM | fT3| , 
cell-cell or cell-obstacle interactions Q and chemotaxis 112|, these models are thus potential 
candidates for the model driving the evolution in our tracking algorithm. 

We work with diffuse interface representations of the cell membrane to make use of the 
mature theory for the optimal control of semilinear PDEs. One attractive aspect of this approach 
is that, as we do not require sharp interface representations of the cell membrane, it may be 
possible therefore to work directly with the raw experimental image data set without any need 
for segmentation. However the diffuse interface or phase field framework we employ does make 
fhe algorifhm compufafionally intensive as evidenced by fhe relafively large CPU limes for our 
experimenls and a key area for fulure work is lo invesligale improvemenls in fhe compulafional 
efficiency of fhe algorifhm. This need is especially evidenf if one wishes fo frack cells in 3d, 
as allhough our Iheorelical framework applies equally lo Ihis selling fhe compulafional cosl 
becomes prohibilive. Compulafional aspecls under invesligalion include 

• Spalial and temporal adaplivily which is challenging in Ihis selling as Ihe solution of Ihe slate 
equation enters Ihe adjoinl equation. 

• Alternative update schemes for Ihe conlrol lo Ihe simple yel robusl gradienl based update 
considered in Ihis sludy. 

• Parallelisation and Ihe developmenl of fasl solvers for Ihe solution of Ihe slate and adjoinl 
equations. 


Our initial numerical investigations suggesl lhal wilh a combination of Ihe techniques oullined 
above il is possible lo efficienlly Irack 3d cell migration and we reporl on Ihis elsewhere. Olher 
potentially attractive directions for fulure work would be lo consider higher order finite elemenl 
spaces for Ihe discretisation of Ihe forward and adjoinl problems or Ihe use of speclral elemenl 
melhods, bolh of which may allow a more accurate solution of Ihe forward and adjoinl problem 
wilh fewer degrees of freedom, hence reducing Ihe memory requiremenls. 

Investigating Ihe performance of Ihe algorifhm wilh real biological dala for differenl cell 
lypes and in differenl environment is an imporlanl and worlhwhile lask. We are currenlly 
applying Ihe algorilhm lo Ihe fracking of in vivo neulrophil migration and intend lo reporl on 
Ihis elsewhere. As mentioned previously one inlerprelalion of Ihe forcing j]* is lhal il account 
for bolh prolrusive forces generated by polymerisation of aclin al Ihe leading edge of Ihe cell 
logelher wilh conlraclile forces generated by Ihe action of myosin motors al Ihe cell rear. Thus 
a potential avenue for assessing Ihe plausibilily of Ihe cell Iracks computed wilh our algorilhm 
would be to compare Ihe computed Vj* wilh experimenlal imaging dala on Ihe location of poly¬ 
merised aclin and myosin-II on Ihe cell membrane wilh Ihe expeclalion being lhal regions in 
which Ihe computed forcing rj* is positive would correspond to regions rich in polymerised 
aclin and regions in which Ihe computed forcing rj* is negative would correspond to regions 
rich in myosin-II. There are also many extensions of our approach which are likely to prove 
useful in applications. Our algorilhm could equally be applied to Ihe identification of (possibly 
lime-dependenl) parameters in models for cell migration (e.g., a spatially conslanl forcing or 
material parameters such as surface tension or bending rigidity) however in Ihis case il is likely 
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that the sharp interface approach we propose in Q will be more efficient. As observed in 
some of the experiments we report on, the framework we employ allows changes in topology 
of the cells. Whilst this may be desirable for some applications, e.g., tracking cells beyond cell 
division or cell fusion, in many biological experiments the topology of the cells is fixed. Our ex¬ 
periments suggest that topological changes arise primarily in the case of multi-cell image data 
sets. In this setting it should be possible to track the evolution of certain topological invariants 
(or more specifically diffuse interface representations of such invariants) and use these as an 
indicator for when the computed cells are changing in topology. The user could then manu¬ 
ally reduce the multi-cell tracking problem to multiple single cell (or smaller scale multi-cell) 
tracking problems by specifying the correspondence between cells in different frames, with the 
hope that changes in topology do not occur for these new problems. The model we propose 
for the evolution in this study is a simplification of more general physically relevant models 
in which bulk or surface PDEs for the biochemistry are coupled to a geometric evolution law 
for the motion. An important area for future work is the extension of the framework to this 
more general setting. We note that the phase field approach we employ makes it computation¬ 
ally straightforward to couple the geometric evolution law for the motion to bulk PDEs (posed 
either within the cell or in the extra-cellular matrix) |[8| [T0||T^ . 

We hope that our optimal control-model fitting based framework is a useful first step to¬ 
wards incorporating advances in the modelling of cell migration into cell tracking algorithms. 
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A Numerical solution of the forward and adjoint prob¬ 
lems 


A.l Discretisation of the state equation 

We introduce the variational form for the forward problem (12.81) defined as follows. Eind 
(v 7,A) e L2([0,r];Hi(D)) xL^{0,T) such that 

f dt(p'iljdx+ f Vip ■ Vijjdx = - f {ccrj — \)il)dx — ^ f G' {ip)'il}dx Vy) G H^(U). 
Jn Jn ^ Jn Jn 

Eet T be a decomposition of U into simplexes (for simplicity we assume Q is polygonal). We 
define the finite element space 


V := e n C0(D) : e Vfe G T}. 


(A.l) 


Eor the time discretisation employ an implicit-explicit method where the diffusive term is 
treated implicitly and the reaction terms explicitly. Introducing the shorthand for a time dis¬ 
crete sequence f"' := f{G) and a uniform timestep r with T = Mr, M G N, the fully discrete 
scheme reads, forn = 0,..., M — 1, given ypjj, G V find , A”"''^) G V x M such that 

- [ - ^h^hdx + f Vip'P^^ • Viphdx = 

T Jn Jn 

^ [ {cGVh - - 4 / ^^{G'{(pl))'iphdx y'lph G V, 

^ Jn ^ Jn 
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where : C'‘’(0) —)• V denotes the Lagrange interpolant. 

We solve the above problem using the iterative technique introduced and studied in |211, 
which uses a bisection method for the Lagrange multiplier. In particular we seek an iterative 
sequence where solves 


- / - ^Vli’hdyi + f ■ Viphd:x. = 

T' Jn Jn 

^ Jn ^ Jn 

with A”+A1 = -^ + 1 , A’"+ 1’2 = ^ _ 1 and {A”+iA+i}fc >2 satisfying 


where we recall 


A/"+‘ - 


(/fibr‘'V - /nbr'-*-'i+) 

;= f |V?J+ + ([»»oiJ+ - b2l+) dx. 


We deem this iteration to have converged when |a"'+^A+i _ ^n+i,fc| < (ol. 

The discretisation of the forward problem without the volume constraint is as above with 
A = 0. The discretisation of the adjoint problem without the volume constraint is the same as 
above as the Lagrange multiplier X{t) does not enter the adjoint problem (2.151. 
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